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Abstract — We derive optimal filters on the sphere in the 
context of detecting compact objects embedded in a stochastic 
background process. The matched filter and the scale adaptive 
filter are derived on the sphere in the most general setting, 
allowing for directional template profiles and filters. The perfor- 
mance and relative merits of the two optimal filters are discussed. 
The application of optimal filter theory on the sphere to the 
detection of compact objects is demonstrated on simulated data. 
A naive detection strategy is adopted, with an initial aim of 
illustrating the application of the new filter theory derived on 
the sphere. Nevertheless, this simple object detection strategy is 
demonstrated to perform well, even at low signal-to-noise ratio. 

Index Terms — Filtering, spheres, convolution. 



I. Introduction 

THE detection of compact objects embedded in a stochas- 
tic background is a problem experienced in many 
branches of physics and signal processing, and as such it has 
received considerable attention. Many of these applications are 
restricted to flat Euclidean space, such as the one-dimensional 
line or the two-dimensional plane. However, data are often 
measured or defined on other manifolds, such as the two- 
sphere. For example, applications where data are defined on 
the sphere are found in planetary science, geophysics, com- 
puter vision, quantum chemistry and astrophysics. Astrophys- 
ical applications include observations made on the celestial 
sphere of the cosmic microwave background (CMB) (e.g. [1]), 
a relic radiation of the Big Bang. CMB observations may be 
contaminated by localised foreground emission due to point 
sources or Sunyaev-Zel'dovich (SZ) effects induced by the 
hot intergalactic gas bound to clusters of galaxies [2]. These 
foreground emissions must be separated from CMB observa- 
tions in order to provide cleaned CMB data for cosmological 
analysis or simply so that they may be studied in their own 
right. Furthermore, other primordial physical phenomena may 
imprint localised signatures in the CMB that are of interest 
(e.g. cosmic strings [3]). The extension of optimal filter theory 
to the sphere would allow compact objects embedded in full- 
sky CMB data to be detected in an analogous manner to that 
performed in the plane currently, i.e. by filtering the observed 
field to enhance the contribution of embedded objects relative 
to the stochastic background, before attempts are made to 
recover the embedded objects from the filtered field using 
various classification schemes. 

When adopting the filtering approach to object detection a 
number of criteria may be imposed so that the filter kernels are 
in some sense optimal. The optimal matched filter (MF) has 
been used extensively in many branches of physics and signal 
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processing, and in the context of detecting point sources and 
SZ emission in CMB observations made on small patches of 
the sky, which are assumed to be approximately flat, by [4] and 
[5] respectively. Other optimal filters in Euclidean space, such 
as the scale adaptive filter (SAF), have been derived by [6], [7]. 
It is shown by [6] that the Mexican hat wavelet on the plane 
is in fact the optimal SAF for the special case of detecting a 
Gaussian source embedded in a white noise background. The 
Mexican hat wavelet has been used to detect objects embedded 
in CMB data on small patches of the sky [8]. Furthermore, 
both the MF and SAF (including the Mexican hat wavelet) 
have been applied to simulated time ordered CMB data to 
detect point sources [9]. Some debate exists, however, over 
the advantage of the SAF relative to the usual MF [10]. 

All of the works discussed previously are limited to Eu- 
clidean space. To analyse full-sky CMB observations the tech- 
niques described previously must be extended to a spherical 
manifold. As a consequence of a full-sky analysis, it should be 
noted that the statistical properties of the background process 
are assumed to be stationary over the sphere. The Mexican hat 
wavelet analysis has been extended to the sphere and applied 
to point source detection in the CMB by [11]. The extension 
of optimal filter theory to the sphere has been derived recently 
by [12] and applied to simulated data to detect the SZ effect 
[13]. However, the optimal filter theory derived by [12] is 
restricted to azimuthally symmetric source profiles and filters 
on the sphere. In this paper we re-derive optimal filter theory 
on the sphere, making the extension to the more general class 
of non-azimuthally symmetric source profiles and filters. In 
addition, we generalise to L p -norm preserving dilations in 
order to highlight some minor amendments to previous works. 

The remainder of this paper is organised as follows. In 
Section HJ some mathematical preliminaries are presented be- 
fore the object detection problem is formalised in Section Hill 
Derivations of the MF and SAF on the sphere are presented 
in Section [IVJ In Section [V] the new optimal filter theory is 
applied to detect objects embedded in simulated data, in the 
MF case. Concluding remarks are made in Section [VD 



II. Mathematical Preliminaries 

It is necessary to outline some mathematical preliminaries 
before attempting to derive optimal filters on the sphere. 
Firstly, harmonic analysis on the two-sphere S 2 and on the 
rotation group SO (3) is reviewed. By making all assumptions 
and definitions explicit we hope to avoid any confusion over 
the conventions adopted. A dilation operator is then defined on 
the sphere, before we review filtering on the sphere. Dilation 
and filtering are described in both real and harmonic space. 
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A. Harmonic representations 

We consider the space of square integrable func- 
tions L 2 (S 2 , dfl(u)) on the unit two-sphere 8 2 , where 
dfl(uj) = smOdOdcj) is the usual rotation invariant mea- 



sure on the sphere and uo 



(0,0) G § 2 de- 



notes spherical coordinates with colatitude and longi- 
tude (j). A square integrable function on the sphere / G 
L 2 (S 2 , dVt{uo)) may be represented by the spherical har- 
monic expansion /M = E^oEl=-J^mH, where 
the spherical harmonic coefficients are given by the usual 
projection on to the spherical harmonic basis functions: 
jtm = J§2 f(v) Yimi 00 ) dfi(o;). The * denotes complex con- 
jugation. We adopt the Condon- Shortley phase convention 
where the normalised spherical harmonics are defined by [14] 



Y tm {u) = (-1)' 



/21+1 (l-ro)! 
47r (£ + m)\ 



PP(cos9)e 



\m(f) 



where P™(x) are the associated Legendre functions. Using 
this normalisation the orthogonality of the spherical harmonic 
functions is given by 



/ 



Y £rn (u;) Yg, m ,{u3) dQ(u) = 6e£>6 n 



(1) 



where Sij is the Kronecker delta function. 

To perform filtering on the sphere one must define trans- 
lations on the sphere, which may be naturally represented by 
rotations. Rotations on the sphere 7Z are characterised by the 
elements of the rotation group SO (3), which we parameterise 
in terms of the three Euler angles p = (a,/?, 7) G SO(3), 
where a G [0,2tt), f3 G [0, tt] and 7 G [0,2tt)Q The 
rotation of / is defined by [R(p)f](co) = f{R~ 1 uj), where 
R p is the rotation matrix corresponding to 7Z(p). It is also 
useful to characterise the rotation of a function on the sphere 
in harmonic space. The harmonic coefficients of a rotated 
function are related to the coefficients of the original function 
by 

t 

mp)fhm= E D e mm ,( P ) f em , . (2) 

m'=-£ 

The Wigner functions may be decomposed as [14] 

D e mm , (a, p, 7) = e~ ima d l mm , (/?) e~ im 'i , (3) 
where the real polar d-matrix is defined by [14] 



, (R , / v J + m')!(l-m')! f . /T m ~ m 

d mm>(P) = * I — — — ( Sin - 



(I + m)\(l - m)\ V 2 

ra'+m 



/ n\ m -\-m 

x^coslj p/r m 7 m ' m+m) (cos/3), (4) 

where P^ a,b \-) are the Jacobi polynomials. The Wigner func- 
tions satisfy the orthogonality condition 



SO(3) 



D f mn (p)D f J n ,(p) dg(p) 



8tt 2 
2t+ 1 



' $mm' $nn' •> (5) 



x We adopt the zyz Euler convention corresponding to the rotation of a 
physical body in a fixed co-ordinate system about the z, y and z axes by 7, 
(3 and a respectively. 



where dg(p) = sin j3 da dj3 d^. Recursion formulae are avail- 
able to compute rapidly the Wigner d-matrices (see e.g. [15]). 

B. Dilation on the sphere 

To perform filtering on the sphere at various scales a 
spherical dilation operator must be defined. We adopt the 
spherical dilation operator first defined by [16] to derive a 
continuous wavelet transform on the sphere. In this case, the 
stereographic projection is used to project the sphere on to the 
plane (see Figure [T]). It is shown by [17] that the stereographic 
projection operator is the unique unitary, radial and conformal 
diffeomorphism between the sphere and the plane. Dilations 
on the sphere are defined by: (a) first lifting the sphere to the 
plane using the stereographic projection; (b) performing the 
usual Euclidean dilation in the plane; (c) reprojecting back on 
to the sphere using the inverse stereographic projection. The 
spherical dilation operator derived by [16], [17] preserves the 
L 2 -norm of functions. We generalise to L p -norm preserving 
dilations in order to highlight the impact of this choice 
on the final optimal filter equations derived in Section HVl 
Although the choice of p is not of considerable practical 
interest, different works have implicitly assumed different p 
values, hence the explicit dependence given here should allow 
a direct comparison with previous works. Definitions of the 
dilation and inverse dilation (contraction) operators follow, 
accompanied by short proofs to show that these operators do 
indeed preserve the L p -norm. 

Definition 1: The L p -norm preserving spherical dilation is 
defined by 

/*,p(M) = [v P {R)f]{0A) = lHR,o)] 1/p f{0i/nA) , (6) 

where scale R G R+, L^-norm p G Z+, tan(0fl/2) = 
i?tan(0/2) and the \(R,6) cocycle term introduced to pre- 
serve the appropriate norm is defined by 

4R 2 

A(jR ' 0) = [(i? 2 -l)cos# + 0R 2 + l)] 2 ' (?) 

□ 

Proof: See Appendix U ■ 
A contraction, or inverse spherical dilation, may be similarly 
defined and follows trivially from Definition 1. 

Definition 2: The L p -norm preserving inverse spherical di- 
lation is defined by 

[V-\R)f)(6^) = [\(R,0 R T 1/p f(0 R ,<l>) . (8) 

□ 

The spherical dilation and contraction operations are per- 
formed about the north pole. A dilation about any other point 
on the sphere may be performed by rotating that point to the 
north pole, dilating, then rotating back to the original position. 

When constructing optimal filters on the sphere it is neces- 
sary to consider the spherical harmonic coefficients of dilated 
functions. We derive an intermediate result here to express 
the harmonic coefficients of the dilated function in terms of 
the original undilated function and dilated spherical harmonic 
functions. 

Lemma 1: The harmonic coefficients of a dilated function 
may be given either by projecting the dilated function on 
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South pole 



Fig. 1. Stereographic projection of the sphere on to the plane. 

to each spherical harmonic, or equivalently by projecting 
the original undilated function on to contracted spherical 
harmonics, scaled by the appropriate cocycle factor: 

Indeed, for the 2-norm case (p = 2) one finds 1 — 1/p = 1/p 
and the spherical harmonics are contracted in accordance with 
the definition of the inverse dilation: 

\W)f] tm = I f(e,4>)[v^(R)Y; m }(e,4>) dfl(6U). 

□ 

Proof: Consider the spherical harmonic coefficients of a 
dilated function 

[ T >p(R)f\i m = I [A(it!^)] 1/ ^(^i/ i? ,0)y; m (^0)do(^0). 

By performing a change of variables this may be represented 

by 

from which © follows by noting dfi(0R,0) = 
Or)} -1 dQ(0 ,(/)). ^ ■ 

C. Filtering on the sphere 

Filtering on the sphere is performed analogously to filtering 
on the plane; it is the spherical convolution of the filter kernel 
with the analysed signal. The analogue of translations on the 
sphere are rotations, thus the filtered field is given by the 
directional spherical convolution 

F(p,R,p)= [ f((j) [K(p)*R,p]*(w) dfi(w) , (10) 

where ^ G L 2 (S 2 , dQ(u)) is the filter kernel. All orientations 
in the rotation group SO (3) are considered, thus directional 
structure is naturally incorporated. The filter equation (171)1) is 
identical to the analysis formula of the continuous wavelet 
transforms derived on the sphere by [16], [17], hence our 
fast algorithm to evaluate this equation [18] may be applied 



to compute the filtered field rapidly. For an example of the 
directional spherical filtering operation applied to Earth data 
see [18]. 

When deriving optimal filters on the sphere it is often 
convenient to represent the filtering operation in harmonic 
space. Representing the analysed function and the rotated filter 
kernel in harmonic space, one may rewrite (ITOb as 

oo £ £ 

F(p,R, P ) = j2 E E hmD"* ml { P )^ R , P ); ml . (id 

£=0 m=-£ m' = -£ 

For a filter centred on the north pole the harmonic represen- 
tation of the filtering operation reduces to 

F(0,R,p)=J2^rn^R, P T im : (12) 

£m 

where here and subsequently we use the shorthand notation 
^2eLoY^m=-e = Ysim- Furthermore, for the special case 
of a filter kernel that is azimuthally symmetric, the filter is 
dependent on 6 only (and not <j)) and, consequently, the filtered 
field is independent of the value of 7. In this case, the spherical 
harmonic coefficients of the filtered field for a particular scale 
are simply given by 

/ 47T 

[F(R,p)] tm = y^rpy hm {*R, P Y m ■ (13) 

III. Problem Formalisation 

In this section we formalise the problem of detecting 
compact objects embedded in a stochastic background noise 
process, and propose optimal filtering to enhance the detection 
of such objects. The formulation given here is similar to that 
of [12] but is considered in the most general sense, allowing 
asymmetric templates of various amplitude, position and orien- 
tation. Removing the assumption of an azimuthally symmetric 
template, the case considered in previous works [12], [13], 
introduces a number of complications as the spherical filtering 
operation may no longer be represented in harmonic space 
simply as a product of spherical harmonic coefficients. 

A. Formulation 

Consider an observed field on the sky f(ou) consisting of a 
number of sources si{uj) embedded in a stochastic background 
process n(u). Observations are likely to be made in the pres- 
ence of additional instrument noise r(u) also. The observed 
field is obtained by measuring (convolving) the actual field 
with some beam b(oj). The beam response and instrumental 
noise may be absorbed into the template and the statistical 
properties of the background, therefore these components may 
be included trivially when required. Without loss of generality, 
the observed field may therefore be represented by 

/H = ^ S ,H+nH. (14) 

i 

Each source may be represented in terms of its ampli- 
tude Ai and source profile Si(u) = Ai Ti(u), where 
Ti{uj) is a dilated and rotated version of the source pro- 
file t(uj) of default dilation centred on the north pole, 
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i.e. Ti(uj) = lZ(pi)V(R ij p) t(uj). One wishes to recover the 
parameters {Ai^Ri^pi} that describe each source ampli- 
tude, scale and position/orientation respectively. The stochas- 
tic background process is assumed to be a zero-mean 
E [n{uj)] = 0, homogeneous and isotropic Gaussian random 
field, fully characterised by the spectrum 

E [n£ m ri£, m t] = S££>5 mrn > Q , (15) 

where E[-] denotes the expectation operator. To facilitate the 
detection of compact sources, the observed field is filtered 
using ([TOt to enhance the source contribution relative to the 
background noise process. The choice of filter kernels that are 
in some sense optimal is addressed next. 

B. Filter constraints 

Various optimal filters may be defined by imposing different 
constraints on the filtered field. Without loss of generality, 
we derive optimal filters for the detection of a single source 
located at the north pole (hence we drop the i subscript that 
denotes source index). Sources located at other positions and 
orientations are found by rotating the optimal filter, i.e. by 
considering filtered field coefficients over SO (3). 

The following filter characteristics may be imposed when 
constructing optimal filters: 

(i) Unbiased: The filtered field is an unbiased estimator 
of the source amplitude at the source position, i.e. 
E[F(0,R,p)]=A. 

(ii) Minimum variance: The filtered field has minimum 
variance at the source position. 

(iii) Local extremum in scale: The expected value of the 
filtered field has an extremum with respect to scale, i.e. 
■^E[F(0,R,p)] =0ztR = R . 

Imposing criteria (i) and (ii) only one obtains the usual 
MF. Imposing the additional constraint (iii) one obtains the 
SAF, first introduced on the plane by [6]. This additional 
constraint imposes an extremum at the unknown scale Rq. In 
the derivation of the SAF this unknown scale drops out and the 
final filter expressions are independent of Ro. The additional 
constraint imposed when constructing the SAF provides a 
means to estimate unknown source scales by checking for a 
maximum in scale but at a cost of reduced gain. This issue is 
discussed in more detail in Section IIV-DI 

C. Filtered field statistics 

In order to derive optimal filters on the sphere it is necessary 
to determine first expressions for the mean and variance of the 
filtered field at the source position. The filtered field mean at 
the source position is given by 

E [F(0, R, p)}=A^ Urn (*fl, P )L • (16) 

ira 

The filtered field variance at the source position is given by 

4(0, R,p) = E [\F(0,R,p)\ 2 ] - \E[F(0,R,p)] | 2 . (17) 
The first term becomes 

E[\F(0,R,p)\ 2 } = \E[F(0,R,p)]\ 2 + J2Ce \{^R, P ) im \ 2 , 

ira 



where we have relied on the fact that the stochastic noise 
process has zero mean and is homogeneous and isotropic. The 
variance is therefore given by 

o 2 F (0,R,p)=J2d \(*R, P ) im \ 2 • (18) 

ira 

The filtered field variance at the source position is also used 
to determine the error on amplitude estimates. 

IV. Optimal Filters 

In this section we derive the MF and SAF on the sphere 
for an arbitrary template profile. The extension of optimal 
filter theory to the sphere has been derived already by [12] 
for the special case of azimuthally symmetric source profiles. 
We re-derive optimal filter theory on the sphere here, making 
the extension to the more general class of non-azimuthally 
symmetric source profiles and filters. In addition, we gen- 
eralise to L p -norm preserving dilations in order to highlight 
some minor amendments to previous works. We show that the 
resultant optimal filters reduce to the expected definitions for 
azimuthally symmetric template profiles and also that, in the 
flat, continuous limit, these forms reduce to the optimal filter 
definitions derived on the plane. To conclude this section we 
compare the performance of the MF and SAF on the sphere 
and discuss the relative merits of each filter. 

A. Matched filter 

Theorem 1: The optimal MF defined on the sphere is ob- 
tained by imposing criteria (i) and (ii) defined in Section HTl-Bl 
i.e. by solving the constrained optimisation problem: 

min o-p(0,R,p) 

w.r.t. (*R, P ) £m 

such that 

E[F(0,R,p)]=A. (19) 

The spherical harmonic coefficients of the resultant MF are 
given by 

= ^ > (20) 

where 

a = Y, C i 1 \ T t™\ 2 • (21) 

ira 

□ 

Proof: See Appendix HH ■ 
A measure of the capability of an optimal filter to detect 
an embedded source is given by the maximum detection level, 
defined by 

r s «M. (22) 

CTF,min(0, 

The minimum variance of the filtered field is found by 
substituting the expression for the MF given by (El) into ([[8]). 
One finds that the MF filtered field minimum variance is given 
by a^ min (0, = a -1 , thus the maximum detection level 
for the MF is 

r MF = a 1 / 2 A . (23) 
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B. Scale adaptive filter 

To construct the optimal SAF denned on the sphere we first 
recast criterion (iii) in a form expressed in terms of the filter 
spherical harmonic coefficients and not the coefficients of the 
differentiated filter. A suitable form is given by the following 
lemma. 

Lemma 2: Optimal filter criterion (iii), imposing a local 
extremum in scale, may be recast in the following more 
applicable form: 



C. Azimuthally symmetric case 

In this subsection we show that the expressions derived 
above for arbitrary template profiles reduce to the forms 
expected for azimuthally symmetric templates (the case con- 
sidered by [12]). The spherical harmonic coefficients of az- 
imuthally symmetric functions on the sphere are non-zero for 
m = only, hence the general filter expressions that we derive 
for a directional template should reduce to the symmetric result 
when setting m = 0. The definition of the filter variables used 
herein differ slightly to those defined by [12] (which we denote 



V (* Ro , p y e JA epTem - B tm r t - lim ) = , (24) b y & '> c ' and A ') ; The relationships between these sets of 



where A ip =l + 2/p - 1 and B im (£ 2 - m 2 ) 1 / 2 . □ 
The proof of this lemma may be found in [19] (the proof is 
not repeated here since it requires multiple pages). It is now 
possible to derive the following theorem, the definition of the 
optimal SAF defined on the sphere. 

Theorem 2: The optimal SAF defined on the sphere is 
obtained by imposing criteria (i), (ii) and (iii) defined in 
Section IIII-Bi i.e. by solving the constrained optimisation 
problem: 

4(0, R,p) 



mm 

w.r.t. (^R Q , p ) i 



such that 



and 



E[F(0,R,p)]=A 



—E[F(0,R,p)} 



(25) 



(26) 



R=Rq 



The spherical harmonic coefficients of the resultant SAF are 
given by 

/ lXr \ CT£ m - b(A£ p T£m - BimTi-i^m) 



where 



Aft 



C = ^ C £ 1 \A£ p T£m — B£mT£- 



l.m 



A = ac-\b\ 2 



(29) 



(30) 



□ 



and a is defined by (LTTT) . 

Proof: See Appendix HIT1 ■ 
The minimum variance of the filtered field is found by 
substituting the expression for the SAF given by (l27l) into (Tl~8T> . 
One finds that the SAF filtered field minimum variance is given 
by crp min (Ro,p, 0) = cA _1 , thus the maximum detection 
level for the SAF is 



SAF 



-l/2 A l/2 Am 



variables are stated in [19]. For now we simply note that all 
filter variables are identical for p = 2, however the case p = 1 
is adopted by [12]. Simplifying the formula for the MF given 
by (l20l) to the azimuthally symmetric setting, one finds that 
the resulting expression is identical to that derived by [12]. 
Simplifying the formula for the SAF given by (l27l) in a similar 
manner, one obtains the following expression for the SAF of 
a symmetric template: 



(* H|P U = (A'Q)- 



T£0 



d\n£ 



It is apparent that the Aq p coefficient of the first b' and of 
a' are in general dependent on the choice of the L p -norm 
preserved by the dilation. For p = 1, the case considered by 
[12], one finds A$\ = 1, not two as given in [12]. When we 
repeat the derivation of the SAF for the azimuthally symmetric 
case given by [12] explicitly, we also get a unity term rather 
than a factor of two for these coefficients, and hence correct 
an algebraic error in [12]. The forms derived herein for the 
MF and SAF on the sphere of a directional template therefore 
reduce to the correct forms derived directly for an azimuthally 
symmetric template. Moreover, in the flat, continuous limit, 
the azimuthally symmetric optimal filters derived on the sphere 
reduce to the forms derived previously on the plane (to make a 
comparison see, respectively, e.g. eqn. (21) of [5] and eqn. (10) 
of [6]). 



D. Comparison 

The relative merits of the MF and SAF are compared in this 
subsection. Recently, the advantages of the SAF (on the plane) 
have been questioned [10], although these concerns have been 
refuted by the original proponents of the SAF [6], [7]. We 
hope to clarify this debate by suggesting that in the ideal case 
the SAF may not provide a theoretical advantage, however in 
practice the SAF may indeed prove advantageous. 

The SAF filter imposes an extremum in scale in the filtered 
field and thus must satisfy an additional constraint to the MF. 
Consequently, the gain of the SAF must be lower than that 
of the MF. One may show this analytically by comparing the 
ratio of the detection levels derived in the preceding sections 
for each optimal filter: 



J- SAF 

Tmf 



= \ 1 



|ft| 2 



ac 



(32) 



(31) 



Filter variables a and c are real and are strictly positive always, 
consequently (f32l) is less than one always. 
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As the gain of the SAF is lower than the gain of the MF 
one would hope that the SAF provides some other advantage. 
This is indeed the case when the scale of the source template 
is unknown. When the source size is unknown the observed 
field may be filtered at a range of candidate scales. It is then 
possible, at the source position, to trace the value of the filtered 
field over scales. The SAF imposes a peak (with respect to 
scale) in the filtered field that may then be related, either 
analytically or numerically, to the unknown size of the source 
template. The drop in gain provided by the SAF is therefore 
offset by the ability to determine the unknown size of the 
source. 

It has been argued that it is also possible to estimate an 
unknown template size using the MF [10], in which case the 
SAF would not provide any advantage. For the MF, it may be 
possible to relate the filtered field curve (with respect to scale) 
to the unknown scale of the template. By fitting the entire 
curve one could therefore directly estimate the source scale, 
although this curve is unlikely to have a local peak. In practice 
it is much easier to determine a peak in the filtered field (the 
SAF approach to estimating source size), than try to fit a curve 
to the field (the MF approach to estimating source size). The 
amplitude of the source is unknown also, hence it is only the 
curvature of the curve that one may fit and curvature changes 
much more rapidly about an extremum. It is therefore clear 
that the SAF does indeed provide some practical advantages 
to the MF when the scale of the source profile is unknown. 

V. Simulations 

We apply the optimal filter theory derived in the preceding 
section to the detection of compact objects from simulated 
data. A very simple detection strategy is adopted. The devel- 
opment of more sophisticated detection strategies, a rigorous 
quantification of the performance of the resulting object de- 
tection and the application to real data is the subject of future 
work. The motivation here is to demonstrate the theory with 
a very simple example only. 

A. Optimal filters 

We construct examples of optimal filters defined on the 
sphere using the filter expressions derived in Section HV] Our 
ultimate aim is to apply filters to detect compact objects 
embedded in CMB data, such as the recent observations 
made by the NASA Wilkinson Microwave Anisotropy Probe 
(WMAP) [1]. Optimal filters are therefore constructed here 
in an analogous simulated setting. The CMB power spectrum 
that best fits the WMAP data is used to model the stochastic 
background process, with a Gaussian beam of full-width-half- 
maximum (FWHM) of 13' applied. Isotropic white noise of 
constant variance 0.05(mK) 2 is added to reflect the noise 
in WMAP observations. We consider the construction in 
this setting of the MF and SAF for the directional butterfly 
template illustrated in Figure [2] (a) (the butterfly template is 
defined by the partial derivative in one direction only of a two- 
dimensional Gaussian on the sphere; see [18] for a definition)]! 

2 The step of the butterfly template may be used to model the line-like 
discontinuity of Kaiser-Stebbins type cosmic string signatures [3]. 



The resultant MF and SAF are illustrated in Figure [2] (b) 
and (c) respectively. Optimal filters are constructed here in 
the context of L 2 -norm preserving dilations, i.e. for p — 2. 
Furthermore, in practice the template profiles are assumed to 
be band-limited at £ max , in which case all expressions and 
summations involving i are computed up to ^ max only. In these 
experiments we choose ^ max = 256 since this is a more than 
adequate band-limit to ensure that all structure of the butterfly 
profile is considered. 

A Fortran 90 package called S2FII0 (§ 2 FILtering) has 
been implemented to compute the equations describing the MF 
and SAF defined by Theorem [T] and Theorem [2] respectively. 
This package makes use of our FastCSWI^l package [18] 
to perform fast filtering on the sphere. Additional numerical 
considerations must be taken into account when implementing 
the filter equations since they often involve dividing the 
template spherical harmonics by the power spectrum of the 
background process, which becomes problematic when taking 
the ratio of two very small numbers. These considerations are 
discussed in detail in [19]. 

B. Object detection 

To demonstrate optimal filter theory applied to object de- 
tection a simple example is considered. Although the optimal 
filters we have derived may be used to detect objects of 
unknown size, in this simple demonstration we consider only 
a butterfly template of a fixed, known size. Since the template 
size is known only the MF on the sphere is applied. The MF is 
applied to simulations of the CMB with artificially embedded 
butterfly sources, in order to recover the positions, orientations 
and amplitudes of the sources. Firstly, the simulation pipeline 
used to construct this data is described. The object detection 
procedure is then described briefly, before the results obtained 
from applying this procedure to simulated data are presented. 

1) Simulation pipeline: To test the application of optimal 
filters on the sphere to object detection, simulated data where 
the ground truth is known is analysed. We have implemented 
a Fortran 90 package called COM^l (COmpact eMBedded 
object simulations) to facilitate such simulations. COMB allows 
the user to embed compact functions on the sphere within 
a stochastic background process. The parameters of the em- 
bedded objects are uniformly sampled from some prescribed 
interval, whereas the background process is specified by its 
power spectrum. Functionality is also included to incorporate 
noise and beams. 

An example of maps simulated using COMB is shown in the 
first four panels of Figure [3J In this example the background 
process is described by the best-fit WMAP CMB power 
spectrum, a beam of 13' FWHM is applied and isotropic 
white noise of variance 0.05(mK) 2 is added (the same setting 
in which the optimal filters discussed in Section IV-AI were 
constructed). A map of simulated butterfly sources is shown 
in Figure [3] (c) and is embedded in the resultant simulated 

3 The S2FIL, FastCSWT and COMB packages are available for down- 
load from http : //www .mrao . cam. ac . uk/~ jdm5 7 / These packages 
require the HEALPix (|http : //healpix . jpl . nasa . gov/] [20] and 
FFTW ( |http : //www, f ftw. org/ I libraries. 
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map shown in Figure [3] (d). In this situation the signal-to- 
noise ratio (SNR) of an embedded object is defined by the 
ratio of the peak amplitude of the template relative to the 
root-mean- squared (RMS) value of the CMB background. For 
the butterfly templates it can be shown that SNR « OA A, 
where A is the amplitude of the object (note that the butterfly 
template definition is normalised, with a peak amplitude of 
0.1 3mK for A = 1.0). In the example illustrated in Figure[3]all 
objects have a constant amplitude of A = 1.0, corresponding 
to SNR = 0.40. Moreover, in this simulation all embedded 
objects have the same orientation (7 = 0°); that is, the 
orientation of the objects are aligned to point to the north 
pole. This is a useful test case since it reduces the complexity 
of the subsequent object detection to that of a symmetric 
template profile, while still adopting the more general optimal 
filter definitions. In the following object detection we also 
consider simulations where the orientations of the objects 
vary, however object orientations are only allowed to lie on 
a discrete, uniformly sampled grid with 7V 7 = 5 samples. By 
restricting the allowable orientations to a grid it is necessary 
only to compute the filtered field for a small number of 
orientations. Obviously in practice object orientations will not 
lie on a grid, however in this case one may simply compute 
the filtered field for a higher orientation resolution (or one may 
use steerable filtersQ to extract the orientation corresponding to 
the most dominant feature, thereby effectively considering all 
continuous orientations; this is the topic of future work). An 
orientational resolution of 7V 7 = 5 is sufficient for the simple 
demonstration presented here. 

2) Detection procedure: For the purpose of this demon- 
stration we perform only a naive detection strategy based on 
thresholding the filtered field. Firstly the initial map is filtered 
using the MF to enhance the contribution due to the embedded 
sources relative to the background. The local maxima of the 
filtered field that lie above a certain threshold are used to define 
detected objects. The filtered field is thresholded at a level 
determined by a constant (N a ) times the standard deviation of 
the map. Typically an N a of 2.5 or 3.0 is used. The magnitude 
and position of each local maximum that remains in the 
thresholded field is used to compute the parameters of detected 
sources. This object detection procedure is implemented in the 
S2FIL package. 

This detection algorithm is applied to the simulated data 

4 Steerable template functions yield steerable optimal filters on the sphere 
since the filter equations do not mix m structure. See [17] for a discussion 
of steerable filters on the sphere. 



TABLE I 

Object detection performance for 7 = 0° 



SNR 




Detections 


Number of sources 






Correct 


False 




0.40 


3.0 


10 





10 


0.40 


2.5 


10 


1 


10 


0.25 


3.0 


6 





10 


0.25 


2.5 


8 


2 


10 


0.20 


3.0 


5 





10 


0.20 


2.5 


7 


6 


10 



illustrated in Figure [51 The filtered field is displayed in 
Figure 0(e) and the objects recovered from the local maxima 
of the filtered field are shown in Figure [3] (f). Notice that 
for this simulation at SNR = 0.40 all embedded objects are 
recovered accurately and no false detections are made. We 
now consider more difficult object detection problems. 

3) Results: The example illustrated in Figure [3] demon- 
strates object detection on the sphere for a relatively easy case. 
By considering more difficult detection situations we examine 
the limits of the simple detection algorithm described above. 
Firstly, only a single fixed orientation is considered, before the 
orientation of embedded objects is allowed to vary. 

For objects of known orientation we examine the effect of 
reducing the SNR of the embedded sources on the number 
of successful and false detections. We also experiment with 
thresholding levels N a = 3.0 and N a = 2.5. The results of 
these tests are given in Table H To ensure minimal false 
detections are made one should threshold at N a = 3.0, 
however to improve the completeness of positive detections 
one may use N a = 2.5, at the expense of a small number of 
false detections. 

We now consider objects with an unknown orientation 
laying on a uniform grid of resolution 7V 7 = 5. Allowing 
orientations to vary introduces an extra degree-of-freedom and 
thus makes the object detection problem more difficult. At 
SNR = 0.40, for thresholding levels N a = 3.0 and N a = 2.5, 
of ten embedded objects we recover six and eight objects 
correctly and make zero and one false detection respectively. 
At SNR = 0.25 we recover five objects correctly and make 
one false detection for N a = 3.0 (N a = 2.5 is not appropriate 
for this low SNR). Finally, we demonstrate in this setting 
the recovery of objects at a range of different amplitudes. 
The actual and recovered object parameters are described in 
Table III and may be observed in Figure HI One sigma errors 
on the amplitude estimates are calculated from the variance 
of the filtered field at the source position (as discussed in 
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(a) CMB 



(b) White noise 



(c) Embedded objects 




(d) Simulated sky with objects embedded (e) Filtered field (f) Recovered objects 

Fig. 3. Embedded object simulation and recovered objects. The CMB background shown in panel (a) is simulated in accordance with the CMB power 
spectrum consistent with current observations, noise of variance 0.05(mK) 2 (panel (b)) is added and a beam of 13 / FWHM is applied. The simulated butterfly 
templates (SNR = 0.40) shown in panel (c) are embedded in the the resultant simulated map shown in panel (d). The MF filtered field of the simulated sky 
is shown in panel (e). The local maxima in the filtered field are detected by thresholding at No- = 3.0 to recover the compact objects depicted in panel (f). 
In this example all embedded objects are detected and no false detections are made. (All maps shown here and subsequently are displayed in units of mK.) 



Section HE]). 

It is important to note that the detection strategies demon- 
strated here are extremely naive and are presented merely to 
demonstrate the new filter framework derived on the sphere. 
A more rigorous approach is to perform more sophisticated 
detection classification schemes, such as the Neyman-Pearson 
test. Moreover, for certain classes of template functions which 
are steerable, one may use steerable filters to extract a single 
orientation over the domain of all continuous orientations. This 
is expected to improve considerably the performance of object 
detection in cases of varying source orientation, to the extent 
that one would expect the performance to match that of cases 
where the source orientation is known. These extensions are 
the focus of future work. 

VI. Concluding Remarks 

We have extended the concept of optimal filtering to a 
spherical manifold. Expressions for the spherical MF and SAF 
have been derived for general non-azimuthally symmetric tem- 
plate objects. For the special case of an azimuthally symmetric 
template we have shown that the general results derived herein 
reduce to the forms derived directly in this setting. Moreover, 
we have also shown that in the flat approximation the optimal 
filters derived on the sphere reduce to the corresponding 
optimal filter definitions defined on the plane. 

The focus of this paper is to derive optimal filter theory 
on the sphere, nevertheless we have also demonstrated the 
application of optimal filters to simple object detection. We 
have generated simulations on the sphere of objects embedded 
in a stochastic background process. Using this simulated data 
we have tested a simple object detection algorithm based 
on thresholding the optimal filtered field. This simple object 
detection strategy has been demonstrated to perform well, even 
at low SNR. 



TABLE II 

Actual and recovered object parameters 



Source 




Amplitude 


Object parameters 

Euler angles 


1 


A 


= 1.48 


(a,/3, 7 ) = (15.0°, 108.3°, 72.0°) 




A 


= 1.66 ± 0.20 


(&,/3,7) = (14.7°, 107.9°, 72.0°) 


2 


A 


= 1.06 


(a,/3,7) = (29.2°, 75.4°, 0.0°) 




A 


= 1.26 ± 0.20 


(d,/3,7) = (30.9°, 75.6°, 0.0°) 
(a,/3, 7 ) = (78.4°, 0.7°, 288.0°) 


3 


A 


= 1.43 




Not recovered 




4 


A 


= 1.19 


(a,/3,7) = (94.4°, 109.1°, 216.0°) 




A 


= 1.34 ±0.20 


(d,/3,7) = (92.6°, 109.3°, 216.0°) 


5 


A 


= 0.81 


(a,/3,7) = (107.8°, 99.6°, 144.0°) 




Not recovered 


6 


A 


= 1.13 


(a,/3,7) = (136.0°, 133.9°, 144.0°) 




Not recovered 




7 


A 


= 0.93 


(a,/3,7) = (172.3°, 82.7°, 144.0°) 




A 


= 1.25 ± 0.20 


(a,/3,7) = (171.9°, 81.2°, 144.0°) 


8 


A 


= 0.95 


(a,/3,7) = (241.2°, 137.2°, 72.0°) 




A 


= 1.02 ± 0.20 


(a,/3,7) = (242.8°, 138.0°, 72.0°) 
(a,/3,7) = (317.7°, 172.3°, 72.0°) 


9 


A 


= 1.00 




A 


= 1.76 ± 0.20 


(d,/3, 7 ) = (287.0°, 180.0°, 72.0°) 
(a,/3,7) = (321.7°, 50.7°, 144.0°) 


10 


A 


= 1.16 




A 


= 1.08 ± 0.20 


(a,/3,7) = (321.4°, 52.5°, 144.0°) 


11 


Not present 






A 


= 0.96 ±0.20 


(d,/3,7) = (47.0°, 67.9°, 72.0°) 



In the future we intend to develop more sophisticated 
detection classification schemes using the optimal filtered 
field. Moreover, for steerable template profiles the use of 
steerable optimal filters is expected to improve considerably 
the performance of object detection in cases of varying source 
orientation. We intend to then apply the resulting object detec- 
tion algorithms to real CMB data to search for cosmic string 
signatures, a theoretically postulated but as yet unobserved 
phenomenon [3]. 
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(a) Simulated sky with objects embedded 



(b) Actual objects 



(c) Recovered objects 



Fig. 4. Embedded object simulation and recovered objects for a range of SNRs, N 1 = 5 and No- = 2.5. See TableUHfor the actual and recovered parameters 
of the embedded sources. 



Appendix I 

Proof of Dilation Definition (Definition 1) 

We prove that the definition of an L p -norm preserving 
spherical dilation does indeed preserve the L p -norm. The 
result for the case p = 2 is given by [16] and an explicit 
proof is given by [17]. For all practical purposes only the 
cases p = {l,2,oo} are considered, nevertheless we prove 
the result for all positive integer p, taking a different ap- 
proach to that of [17]. First consider the case of positive 
integer p < 00. One requires \\V p (R)f\\ p = \\f\\ p , or 
equivalently I\ = I2, where I\ = J g2 |/(0, (j))\ p dfi(0, (j>) 
and I 2 = f g2 I [\(R, f(9i/ R , (j))\ p dQ(6, </>). By making 
a change of variables I2 may be rewritten as 



and 



-/ 



\f(9^)\v\\(R,9 R )\sm9 R R 



COS' 



\0 R /2) 



COS' 



\9/2) 



d(9c 



Thus, for 1 1 = I 2 it is apparent, after a little algebra, that 
the cocycle required to preserve the L p -norm (for positive 
integer p < 00) is of the form given by ©. Now consider the 
L 00 -norm case, corresponding to the case where no cocycle 
term is applied. When no cocycle is applied the function is 
dilated only and is not scaled, hence the maximum absolute 
value of the function over its domain remains unaltered. 
Consequently, the infinity norm is also preserved. 

Appendix II 
Proof of MF Theorem (Theorem 1) 

We solve the MF optimisation problem by minimising the 
Lagrangian 



C = J2 Ci |(^, p ), m | 2 +/iiRe 



+ /12 Ittl 



tm 



where /ii and /i2 are Lagrange multipliers. Notice that the real 
and imaginary parts of the constraint are made explicit. To 
solve this problem the filter and template spherical harmonic 
coefficients are represented in terms of their real and imaginary 
parts: let (V RiP ) £m a^ m + \bt m and n m C£ m + k^ m . To 
minimise the Lagrangian one requires 



dC 

dd£ m 



= 2Ci ae m + /ilQ m + \l2dtm = 



dC 

db £m 



plus the original constraints specified by dT9b (one for the 
real and one for the imaginary component). Solving these 
equations simultaneously, one finds fi\ = — 2a _1 and /i2 = 
for the Lagrange multipliers and ae m = a~ 1 Cl 1 C£ m and 
btm = a~ x C^ x dum for the real and imaginary parts of the 
filter spherical harmonic coefficients, where a is defined by 
(|2T1) . Thus, the spherical harmonic coefficients of the optimal 
MF on the sphere are given by {^ RlP ) irn = a^C^Um. 

The extremum found is checked to ensure it is a minimum. 
The second derivatives of the Lagrangian are C aa = 2Q, 
Cbb = 2Ci and C a b — 0, where the subscript notation for 
partial derivatives is used. Now C 2 ah < C aa ^bb and both 
£aa > and Cbb > 0, hence a minimum has indeed been 
found. 

Appendix III 
Proof of SAF Theorem (Theorem 2) 

We solve the SAF optimisation problem by minimising the 
Lagrangian 



tm 



tm 



/i 2 Im (VRo,p)* ern 

.tm 

/i 3 Re ^2 Ro&TimiAipTim - BimTe-i^m) 



/i4 Im 



. tm 



R ,p)* im {Ai p T£ m - BzmTg-^m) 



where /ii, ^2, ^3 and ^4 are Lagrange multipliers. No- 
tice that the real and imaginary parts of each constraint 
are again made explicit. To solve this problem the fil- 
ter and template spherical harmonic coefficients are rep- 
resented in terms of their real and imaginary parts: let 
(^Ro, P )i m &tm + ibe m and r^ m = c irn + id im . To minimise 
the Lagrangian one requires 

dC 



daer, 



n 

■ IJ>4(A£pd£ r 



BtmC£-l,m) 
■ Bimdi-i^m) — 



(33) 
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and 



dC 

db £m 



=2C(bi rn + fiit 



(34) 

plus the original constraints specified by (f25l) and ([26b . Using 
(f33l) and (f34l) it is possible to eliminate a^ m and bi m from the 
original constraints, yielding a system of linear equations for 
the Lagrange multipliers that depend on the template spherical 
harmonic coefficients only. Solving this system one finds 
[/ii,/i2,M3,/M] = f [-c,0,Re(&),-Im(&)] for the Lagrange 
multipliers. Substituting the Lagrange multipliers into (l33t and 
(l34l) one obtains the following expressions for the real and 
imaginary parts of the filter spherical harmonic coefficients: 

aim = (AQ)" 1 [c C£ m - Re(b)(AipC£m - Be m ce-i im ) 

+ Im(6) (A£pd£ m - ^m*-l,m)] 

and 

bim = (AQ) _1 [cd£ m - Re(b)(A£ p d£ m — B^dg-x^m) 
- Im(b)(Ai p ci m - B £ 

Thus, combining the real and imaginary parts, the spherical 
harmonic coefficients of the optimal SAF are given by ([27b . 

The extremum found is checked to ensure it is a minimum. 
The second derivatives of the Lagrangian are C aa = 2Ci, 
Cbb = 2Ci and C a b = 0, where the subscript notation for 
partial derivatives is used. Now C 2 ah < C aa £>bb and both 
£ aa > and £55 > 0, hence a minimum has indeed been 
found. 
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